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INTRODUCTION 


This report outlines a computerized approach for the structural analysis of the 
time-independent cyclic plasticity response and of the metal fatigue failure process. 
The approach combines three main analytical components, as follows: 

(a) A cyclic plasticity model which relates the material's uniaxial stress-strain 
behavior to the multiaxial response of any structural component . 

(b) Damage accumulation criteria which indicate both the life to crack initiation 
and the rate of crack growth , up to complete failure , for metallic structural compo- 
nents that undergo local cyclic plasticity strains . The required test parameters are 
derived from only the fatigue life of smooth material specimens when subjected to 
constant uniaxial plastic strain cycles. 

(c) A finite element model for the numerical solution of the structure's nonlinear 
static and dynamic equilibrium equations. The isoparametric finite elements of the 
plane-stress, plane-strain, and axisymmetric types are incorporated. These elements 
are adequate for the representation of the behavior of most aircraft structural com- 
ponents that undergo meaningful plasticity strains . 

The present combined approach enables the following types of analysis: 

(a) The analysis of cyclic plasticity time-independent and rate-independent 
structural response under any varying loading which induces either proportional or 
nonproportional stress variations. Basically, the analysis is related to the material's 
cyclic steady-state behavior; however, the material's cyclic transient behavior can 



also be approximated . The effect of the cyclic yield stress change is not included , 
and the material is assumed to be of the so-called Masing type, which characterizes 
the metallic alloys used in aircraft . In addition , the material is assumed to be 
initially isotropic. The effect of the material's cyclic anisotropy due to the 
Bauschinger phenomenon is incorporated . 

(b) Crack initiation prediction under varying loadings. The prediction is made 
by employing the Coffin-Manson criterion for the multiaxial stress state. 

(c) Crack growth rate prediction. This prediction is made by employing a novel 
damage criterion which relates crack growth rate to the inverse damage gradient 
along the crack path. The criterion accounts for (1) the effects of plasticity, (2) the 
effects of residual stresses and of multiaxial stress redistributions at the crack tip 
which lead to crack retardation, (3) the effects of multiple overloads and negative 
loads, and (4) the interaction of close cracks. The effect of possible crack closure 

is not directly incorporated; however, this phenomenon is approximated by including 
the effect of the residual compressive stresses at the crack tip , which is the main 
cause of crack closure. The effects of loading frequency, temperature, and other 
time-dependent phenomena are not incorporated . 

(d) Propagated crack growth rate prediction. This prediction is based on the 
application of the above-mentioned damage criterion using developed damage data 
accumulated from several updated finite element models . No procedure for the 
inclusion of residual stresses in the propagated crack's wake is included. It is 
assumed that the effect of these residual stresses in the crack's wake is negligible 
because of their usually small magnitude and because of their accelerating relaxation 
rate . The orientation of the propagated crack is set either normal to the computed 
principal tensile stress or in a direction selected by the user upon consideration of 
the direction of the most damaged paths . 

The computer program is an extension of the NONSAP program (ref. 1) . It incor- 
porates cyclic plasticity models and damage accumulation criteria and has an option 
for sorted output. A full listing of the program's new features is given in the 
appendix of this report. The two-dimensional isoparametric finite elements and the 
numerical solution procedures are those of the NONSAP program. As this program 
is an in-core solver, the size of the finite element model is limited. However, the 
analysis of structural components is still practical by using 130K of the computer 
core, as demonstrated later in this report. 


APPLICATION TO DAMAGE -TOLERANT AIRCRAFT DESIGNS 


The damage tolerance requirements specified in MIL-A-83444(USAF) (ref. 2) are 
based on the assumption that a crack already exists in each element of a new structure 
as a result of flaws in the material , corrosion , or manufacturing damage . The struc- 
ture should sustain the growth of these assumed cracks without a total failure during 
its lifetime , and also still sustain a specific residual static strength . Reference 2 
defines two approaches for the substantiation of a structure's damage-tolerant integ- 
rity: the fail-safe approach and the slow-crack-growth approach. The fail-safe 
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approach assumes a smaller initial crack length and a shorter loading spectrum than 
the slow-crack-growth approach; however, it requires more structural accessibility 
for inspection and overall, as well as local, structural redundancies that are fre- 
quently impractical in aircraft structures. The slow-crack-growth approach can be 
applied to all structural types, and it is also simpler to implement. 

The slow-crack-growth approach requires that crack growth be slow enough not 
to achieve an unstable size during the life of the structure . The initial crack length 
is assumed to be on the order of 0.25 inch (6.3 millimeters) (ref. 2) , and the crack 
is also often assumed to be through the member's thickness. These requirements 
can lead , under the usual applied loading , to significant plastic strains at the crack 
tip . A typical loading specturm is composed of varying tension-compression 
components , with multiple overloads , as depicted in figure 1 . Although the loading 
variation is not of a fully cyclic type , it still often imposes cyclic plasticity stresses , 
because of the material Bauschinger phenomenon . 

The imposed cyclic plasticity at the crack tip and the resultant residual stresses 
exclude the implementation of the usual analytical methods , which are based on the 
stress-intensity range. The present computer program can handle these phenomena 
analytically, by combining the finite element method, the material's cyclic plasticity 
model , and the damage accumulation criterion . This analysis is essential both for 
ensuring the integrity of the structural components during their life and for the 
proper evaluation of the results of structural proof tests . 

The present computer program can also be applied to cases in which the crack 
tip undergoes relatively small cyclic plasticity strains. This application can be 
carried out by idealizing the material's stress-strain uniaxial curve with both a low 
yield stress and a first segment's slope which differs only slightly from the material's 
Young's modulus. However, it should be noted that the accuracy of the present 
damage criterion decreases for smaller plastic strains, while the accuracy of the 
simpler stress-intensity range approach increases. The present damage criterion is 
suitable only for cyclic plasticity strains; therefore, monotonically increased plastic 
strains as exhibited in the static residual-strength analysis cannot be handled by the 
present computer program . In addition , repeated loads which do not cause reverse 
plasticity , but cause plastic reloading at the same unloading stress , are assumed to 
contribute to the cumulative plastic strain but not directly to the cumulative damage . 
This will be clarified later in this report. 

The present computer program does not account for the beneficial effects of 
initial compressive stresses due to shot peening, fastener interference, cold-working, 
and the like . However , it should be realized that these effects are usually small 
because of the quick stress relaxation in the cyclic plasticity field. 

The required input data for the computer program are outlined later in this report. 



THEORETICAL APPROACH 


Cyclic Plasticity Models 

Three plasticity models are incorporated in the present computer program . They 
differ from each other in their definitions of the incremental translation of the yield 
surfaces during the hardening of the material. The three models are identical for the 
proportional stress state, but they lead to somewhat different results for the usual 
nonproportional stress state . As none of these models has yet been shown through 
solid experimental evidence to be superior to the others, the choice of model is left 
to the user. 

The three plasticity models are based on classical incremental time-independent 
and rate-independent plastic flow theory for initially isotropic materials. Incremental 
plastic flow theory assumes that the plastic strain increment is much higher than the 
adjacent elastic strain increment, and that plastic strain increments can be computed 
independently on the basis of the previous loading step stresses . Therefore , small 
loading step sizes, specified by the user, are mandatory for solution accuracy. The 
material's uniaxial stress-strain curve can be idealized by a maximum of three 
elastoplastic piecewise linear segments in addition to the first elastic segment, as 
shown in figure 2(a) . The reversal uniaxial segments are shifted by the program to 
twice the initial yield stress , and the length of the segments is magnified by a factor 
of two , assuming material of the Masing type . However , the user can change the 
idealization of the first reversal in order to represent the material's transient 
condition . 

Each linear segment of the material's uniaxial curve is related to a yield surface 
in the multiaxial stress state, as shown in figure 2(b). Each yield surface is defined 
by the von Mises criterion and the associated plastic flow normality rule . It is allowed 
to translate in the stress space up to its bounding yield surface , to which it remains 
connected until the unloading stage . The translation rate is governed by one of the 
following three hardening rules (fig. 3). Prager's hardening rule physically 
assumes that the incremental translation is in the direction of the plastic strain incre- 
ment, i.e. normal to the yield surface. In order to satisfy this rule unconditionally, 
the surfaces' translations in the zero stress directions are mathematically permitted. 
Ziegler's hardening rule assumes that the incremental translation is in the direction 
of the vector which connects the center point of the current yield surface to the 
existing stress point. Both of these hardening rules require continuous position 
corrections of the yield surfaces to ensure tangency among the surfaces in contact. 
Mroz's hardening rule is based on the inherent fulfillment of this tangency require- 
ment. 

The full mathematical expressions of the plasticity models are presented in refer- 
ence 3 . 

It should be noted that the cyclic plasticity room temperature stress relaxation 
phenomenon is not included in the present plasticity models; however, this phe- 
nomenon is directly included in the present damage criteria , which are discussed 
next . 
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Life to Crack Initiation 


According to the presently used criterion , crack initiation occurs after 2N 
reversals of cyclic loading, when the cumulative damage D equals a unit. The damage 
is expressed mathematically as follows (ref. 4): 


D = 



( 1 ) 


The quantity denotes the integration of the equivalent plastic strain incre- 


ment , de^ , through each pair of reversals . The equivalent plastic strain increment 
is a positive scalar composed of the multiplication of the plastic strain increments 


deP. (i , j = 1 , 2 , 
follows: 


3 in tensor notation) , and it is computed by the plasticity model as 


de* 3 = (^deP. deP.V^ 
V i] i]/ 


( 2 ) 


The quantity o m is the average value of the mean stresses at the two plastic un- 
loadings which define the specific pair of reversals , or 


m 


1/2 [( a ii/3)Fi rs t unloading + (°ii^)second unloading] 


(3) 


The quantity o m also represents the effects of the tensile versus compressive stresses. 

If the reversal loading results in a symmetric stress variation, or a.. = a.. . , 

° n, max li, min 

then o = 0 . 
m 


If the stress relaxation effect is to be included , as it should be when a is not 

m 

small, the user must define an experimental material parameter r (ref. 3) such that 

the relaxed a value becomes 
m 


o = o' (2N) 
m m 


-r(y'de^ > /2) 


(4) 


where is the original average mean stress. For the numerical examples to be 

shown later in this report, a value of r of 277 has been adopted for aluminum alloy 
7075-T6 plate. 

The material parameters n' , c, e£, and o^ in equation (1) are defined by the user 
for the specific material. The parameter n' is the material's uniaxial cyclic exponent. 
It relates the uniaxial stress amplitude, Ao/2, to the applied constant plastic strain 
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amplitude, Ae^/2, in the form of Acr/2 = K' (Ae^/2) n , where K' is assumed to be 

n 1 

approximated by a f / (ep . The value of the exponent n' can be derived from several 

uniaxial plastic strain tests at the material’s cyclic steady state, as indicated by 
figure 4(a). The parameter ej. is the material's cyclic ductility parameter, which is 

smaller than the monotonic ductility parameter, e^. The parameter is the material's 

fracture strength. The parameter c is the Coffin -Manson exponent (fig. 4(b)), which 
is derived from constant plastic strain amplitude tests of the material's uniaxial un- 
notched specimens. 

The values of these material parameters depend on the specimen's surface treat- 
ment and environmental conditions . Therefore , the above-mentioned uniaxial tests 
have to be conducted under the same conditions as exist in the real structure. 


Crack Growth Rate 


The crack growth rate is approximated by the inverse damage gradient along the 
crack path. The cumulative damage is computed by equation (1) at two discrete 
points in front of the crack tip . These discrete points are defined by the two inte- 
gration points of the finite element adjacent to the crack tip. Figure 5 designates 
these integration points as number 1 and number 2; they are located at distances 
of a^ and a 2 from the crack tip, respectively. Assume that the accumulated damage 

at points 1 and 2 is termed and D2, respectively. If the crack propagates by the 
small distance of (a„ - a.) , the damage at point 2 becomes D • thus, the average 
cumulative damage value is 1/2 (D^ + Dg). T ^ ie crack growth rate, ) > is approx- 
imated as follows (ref. 3): 


da 

d (224 ) 



(5) 


where a is half the length of the existing crack. Equation (5) indicates that a complete 
fracture occurs when D 2 > D^. 


The finite element integration points , whose cumulative damage values are used 
for the crack growth rate prediction, are chosen by the user according to the pre- 
dicted crack path , which is usually normal to the direction of the principal tensile 
stress. These integration points should be well within the material's cyclic plasticity 
range. This requires a reasonably small finite element to be used at the crack tip. 


Damage Accumulation Technique 

The damage criterion in equation (1) is applied to each pair of reversals sepa- 
rately, and the results are accumulated during the entire applied loading history. 
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Each pair of reversals is defined , as mentioned before , during two subsequent plastic 
unloadings made in reversal directions . The plastic unloadings in figure 6 , for 
example, occur at points B, D, F, II, J, and L. However, the unloading at point F is 
not considered because the following plastic unloading, at point II, is not in the 
reversal direction. Therefore, the first pair-reversal is AB-CD, the second pair- 
reversal is EH-IJ, and so on. 

For tensile loads , the present pair-reversal damage accumulation technique could 
lead to somewhat more conservative results than the well-known rainflow technique 
(ref. 2) . This is because the rainflow technique refers only to closed loops; in 
figure 6 , the plastic strains along the AB , E'F, G'H branches would not be considered, 
because no closing counterpart branches exist. However, the rainflow technique 
does consider the effect of the elastic loop FGG' . 

The present damage accumulation technique does not account for the effect of 
elastic reversals, i.e. it ignores the effect of the elastic loop FGG' in figure 6. This 
is justified because the damage criterion (eq. (1)) employs the material's cyclic 
ductility strain ej., which is smaller than the material's monotonic ductility strain. 

The technique does incorporate the cyclic parameters n' and c; thus, it is assumed 
that the fatigue damage is due mainly to the plasticity cycles . 


Finite Element Modeling and Equation Solutions 

The NONSAP program's two-dimensional isoparametric elements and its solution 
procedures (ref. 1) are utilized. The eight-node element with undistorted shape and 
3X3 integration points has been found to furnish a suitable representation of both 
the plastic strain variation and the damage gradient. The finite element adjacent to 
the crack tip should be small enough for the two integration points along the predicted 
crack path to be well within the cyclic plasticity range . In addition , the idealization 
should be such that the existing crack front is at the corner node , not at the mid- 
node , of the eight node element . Far from the crack tip and far from the stress 
concentration zones , the number of nodes can be reduced to four to save computer 
core and time . 

The behavior of large plastic strains is approximated by employing the Green- 
Lagrange strain tensor and the second Piola-Kirehhoff stress tensor in Lagrangian 
coordinates. The use of this approximation is justified, since most of the fatigue 
failures are accompanied by only small to moderate cyclic strains around the 
material's yield strain. 

The nonlinear equilibrium equations due to the plasticity and the large strains 
are solved incrementally . The size of the loading steps is variable and is set by the 
user, based on his numerical experience. Usually, several short trial and error 
runs are expected for each specific case before the largest possible step sizes are 
determined. The parameter which usually governs the step sizes is the material's 
uniaxial stress-strain slope. A smaller material slope requires a smaller step size. 

The static analysis requires the construction of a new tangent stiffness matrix at 
each loading step . The dynamic analysis can be carried out either by employing 
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Newmark's implicit time-integration method or by employing the explicit central- 
difference method. The central-difference method is much less time consuming, but 
it is more prone to numerical instabilities and thus requires smaller step sizes. This 
method is especially attractive for cases of small material hardening, where the 
required time step sizes are already relatively small because of the small material 
slopes . The iterative NONSAP procedure for equilibrium corrections is not incor- 
porated because of the possibility of nonconvergence at the plastic unloading steps. 


PROGRAM OUTLINE 


The program utilizes the NONSAP computer program's elements and solution 
techniques for large strains and plasticity, and for static or dynamic analysis. The 
new features presented here include the following: 

(a) The incorporation of the cyclic plasticity models and fatigue data computations 
through a separate overlay (number 3.8; see appendix) . The NONSAP overlay tree 

is shown in reference 1 . 

(b) Sorted output data. This is necessary because of the enormous available 
output data and the need to segregate the fatigue data required for the computation 
of the damage criteria . 

Following is a brief summary of the main computation steps. 

• (a) The overall linear stiffness and mass matrices are constructed first . If 
dynamic analysis is required and Newmark's direct time integration technique is 
used . the overall linear effective stiffness matrix is constructed . In addition , the 
applied load vector is constructed. The large strain stiffnesses derived by using the 
Total Lagrangian procedure and the cyclic plasticity stiffnesses are updated at each 
loading or time step. These stiffness values are added to the linear stiffness matrix. 

(b) The equilibrium equations are solved incrementally, and displacements and 
strains are obtained for each step . The program has an optional two-step restart 
capability, which is useful for problems which involve only partially different loads 
and for dividing a long computer run into two separate and more manageable runs. 

(c) For each finite element integration point which is pre-defined by the user as 
an elastoplastic element, the following steps are executed at each loading or time 
step . 

•The previous step's values of elastoplastic stiffness are recomputed. 

• The plastic strain increment is computed, as is the total equivalent plastic strain. 

• The stress increment is computed and the total stresses are updated. The mean 
stress is computed. 

• The yield surface translations are computed, ensuring that the surfaces' non- 
intersection requirement is met. 
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•The elastoplastic stiffnesses are updated in four subincrements and added to the 
overall structural stiffnesses for the next loading step. 

• Continuous checks are made for plastic unloading. If it occurs, the peak 
von Mises stress is kept in the memory to indicate the following reloading state. 

•Plastic loading or reloading is considered when the current stress point reaches 
the first yield surface . The plastic reloading criterion distinguishes between re- 
yielding at the reversed plastic region and reyielding at the same plastic region. Re- 
yielding at the same plastic region is initiated when the accumulated elastic work 
during the unloading range is zero, or nearly zero. The computed accumulated 
damage value is for each pair of reversals; only fully reversed stress cycles are 
considered . 

I 

•The fatigue data for equation (1) are computed. After each pair of reversals, 
damage is accumulated for an indication of the life to crack initiation . The crack 
growth rate is computed by substituting the results of equation (1) into equation (5). 

The mathematical formulations are presented in reference 3 . 


INPUT AND OUTPUT DATA 


The input data are identical to the NONSAP specifications , with the following 
exceptions . The specified material model number for the cyclic plasticity analysis is 
NPAR(15) = 9. The number of constants per property set should be specified as 
NPAR(17) = 15, and the dimension of the storage array should be specified as 
NPAR(18) = 27. Then the material properties are specified on two input cards. The 
first input card contains eight parameters, in 8F10.0 format, as follows: the Young's 
modulus , the Poisson ratio , the yield stress , and the uniaxial slope of the first elas- 
toplastic piecewise linear segment; the yield stress and the uniaxial slope of the 
second segment; and the yield stress and the uniaxial slope of the third segment. The 
second input card contains seven parameters, in 7F10.0 format, as follows: the yield 
stress and the uniaxial slope of the first, second, and third plastic reversal segments; 
and a seventh parameter , RULE , that indicates the required cyclic plasticity model . 

If RULE = 0 , rigid plastic material is assumed . If RULE = 1 , the well-known isotropic 
hardening rule is employed . If RULE = 2 , 3 , or 4 , the kinematic hardening rule due 
to Prager, Ziegler, or Mroz is used, respectively. 

For a material in the cyclic steady state , the specified reversed yield stresses 
and slopes should be identical to the values of the first reversal. Different slopes 
can be specified for the first and second reversals for representation of the material's 
transient state . In the following reversals the data specified for the second reversal 
are used . 

The output data are printed on four tapes: TAPE6, TAPE12, TAPE13, and TAPE14. 
TAPE6 includes the input data and deflections. TAPE12 includes parameters for 
fatigue analysis. Included are the following terms: 



NEL - The finite element number. 

IPT - The integration point number. 

LO - The number of plastic reversals. For the first plastic range LO = 1 , for the 
second plastic reversal LO = 2 , and so on. 

IPEL - The current position of the equivalent von Mises stress. If IPEL = 1, 2, or 
3. the stress point is on the first, second, or third piecewise linear segment, 
respectively . 

DEPC - The cumulative equivalent plastic strain. 

SAIEAN - The mean stress . 

FT - The equivalent von Mises stress. 

SX - The maximum principal stress. 

SY - The minimum principal stress. 

ALPHA - The direction of the maximum principal stress relative to the element's 
coordinates. 

DWE - Numerical stability indicator . It equals the stress increment times the 
elastic strain increment. The value should be positive; otherwise it indicates 
that a numerical instability due to too high step size has been introduced . 

HP - Numerical stability indicator . It should be equal to the input slope of the 
specific material segment. 

WP - Unloading indicator. If WP is negative, unloading occurs. 

IRE - Reloading indicator. If IRE = 0, there is no reloading. If IRE = 1 or 
IRE = 3 , fully reversal plastic reloading occurs . If IRE = 2 , plastic reloading 
occurs at the same unloading point. 

WP2 - The cumulative plastic work. Used for reference. 

DEE - The current total work. Used for reference. 

TAPE13 includes the computed stresses. TAPE14 includes the computed strains, 
surface translations , and other parameters explained in the printout shown in the 
appendix of this report. 

The output data from TAPE 12 are used for the fatigue analysis. The other data 
used in the fatigue analysis include the material's cyclic stress-plastic strain exponent 
n' and the Coffin-AIanson material parameters c, e£, which are defined in 

equation (1). Also needed is the material stress-relaxation exponent r, which is 

defined in equation (4) . The f de^ 1 value in equation (1) is calculated by subtracting 
the computed DEPC values at the two plastic unloading points which define the specific 
pair of reversals. The average of the SAIEAN values at these two unloading points is 
calculated according to equation (3) . This value should be iteratively reduced by 
employing equation (4) because of the assumed cyclic plasticity stress relaxation. 

Then equation (1) is employed for the accumulation of the pair-reversal damage. 

When it reaches a unit value , crack initiation is assumed. The crack growth rate is 
approximated using equation (5) by substituting the cumulative damage values at the 
two discrete points in front of the crack tip and along the predicted crack path. The 
crack growth path is usually predicted to be normal to the direction of the principal 
tensile stress, which is indicated by the ALPHA value. 
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APPLICATION EXAMPLES 


This section describes the application of the present approach to the analysis of 
two structural components: a cracked panel under variable uniaxial loadings and 
stiffened aircraft skin panel under compressive loading. 

The cracked panel is shown in figure 7(a) . The magnitude of the applied loadings 
is such that significant plastic strains develop in front of the crack tip. Figure 7(b) 
depicts the finite element model, which employs plane-stress four-to-eight node 
isoparametric elements . The eighth node elements are solved by 3 X 3 integration 
points . The uniaxial cyclic material curve , idealized by three piecewise linear seg- 
ments, is shown in figure 8. The material's fatigue properties are based on the 
constant strain amplitude test data from reference 5 . The fatigue ductility parameter , 
e£, is assumed to be 0.18, while the measured monotonic ductility, ep is 0.41. The 

fatigue strength, crL is assumed to be equal to the monotonic fracture strength, cr f , or 

2 1 1 
75.9 kg/mm (108.0 ksi) . The Coffin-Manson exponent c in equation (1) is estimated 

to be 0.52. The material uniaxial cyclic exponent , n' , is 0 .11. 

In order to account for the stress relaxation, a value of r of 277 is assumed in 
equation (4). This value causes the relaxation of the existing mean stress down to 
0.01 percent of its initially computed value, within two fully reversed strain cycles of 
O.leJ.. No experimental evidence exists for this value. 

Results for fully cyclic loading and for tensile cyclic loading are shown in 
figures 9(a) and 9(b) and compared to test results which induce only small plasticity. 
These comparisons illustrate the significant crack growth retardation due to the 
plasticity stress redistributions and due to the residual compressive stresses 
developed after plastic unloading. The relative crack growth retardation is more 
significant for the tensile cyclic loading (fig. 9(b)) than for the fully cyclic loading 
(fig. 9(a)) . This is because the residual compressive stresses in the latter case are 
followed by residual tensile stresses which diminish their beneficial effects. The 
computed crack displacements indicate that no crack closure occurs for the present 

loading conditions. The crack growth rate, , in figures 9(a) and 9(b) is 

depicted as a function of the stress intensity range AK = p • Ao n • vs:, where p is a geo- 
metric parameter, Aa^ is the net section stress range, and a. is the half crack length. 

For cases of small and localized plasticity , the stress-intensity range is generally a 
representative parameter. However, in cases of gross plasticity, as in the present 
examples, AK loses its general validity; thus, the results shown in figures 9(a) 
and 9(b) are specific for the crack length used. 

Figure 10 shows the effect of a tensile overload on the crack growth rate as 
computed by the present approach. It is apparent that this effect becomes more 
significant with increasing values of overload . This is in general agreement with 
the test data that have been reported in the literature . 
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The stiffened skin panel is shown in figure 11. The integral stiffeners' cross 
section at the spar location is changed as shown in figure 11 (b) . Axial loads due to 
overall wing bending could lead to high stress concentrations at the indicated point . 
These stress concentrations can usually be significantly reduced by the addition of a 
small area of structural reinforcement. Two cases, with different reinforcement area 
sizes, are analyzed. They are designated case 1 and case 2. Figure 12(a) shows the 
finite element model used. The applied loads are compression and vary with the 
stiffener's depth, as shown. The applied loading variation, shown in figure 12(b) , 
causes local compressive yielding and high residual tensile stresses after unloading. 
Thus, although no tensile loads are applied, a cyclic compression-tension stress-strain 
field exists, causing crack initiation and propagation. The material's uniaxial stress- 
strain curve is idealized by three linear segments, as shown in figure 12(c) . The 
material’s fatigue properties are the same as those indicated for the cracked panel in 
the previous example . 

Figure 13(a) shows the computed damage curves. Each curve indicates the equal 
damage accumulation value. As depicted, the small reinforcement area in case 2 
significantly improves the life to crack initiation. Figure 13(b) shows the von Mises 
equivalent stress distribution for case 1 . It is apparent that the stress gradient is 
much smoother than the damage gradient. This demonstrates the inability of stresses 
to predict the fatigue failure in a plastic field . 

Figure 14(a) shows examples of the used cracked finite element models. The left- 
hand model represents the initial crack pattern, which is perpendicular to the com- 
ponent's free edge (and to the direction of the principal tensile stress) . However, in 
order to maintain the element's parallelogram shape, which is an important factor for 
numerical accuracy, the crack's direction is changed slightly, as shown. The right- 
hand model in figure 14(a) represents progressive crack growth. The damage curves 
before the crack changes its direction are shown in figure 14(b). The damage 
accumulation gradient and the crack growth rate are derived from the curves shown 
in figures 13(a) and 14(b). The results are summarized in figure 14(c). 


CONCLUDING REMARKS 


This paper describes a computerized approach to the calculation of cyclic plasticity 
structural response , the prediction of life to crack initiation , and the prediction of 
crack growth rate. The method uses three analytical items: the finite element method 
and its associated numerical techniques for nonlinear static and dynamic analysis , the 
material cyclic plasticity theory, and the cumulative damage criteria. 

The required input data include the loading spectrum, the material's cyclic 
uniaxial stress-strain curve, the material's cyclic stress-plastic strain exponent, and 
the Coffin-Manson low-cycle fatigue parameters. These parameters are derived from 
only smooth uniaxial specimens. The method also requires the material's stress 
relaxation exponent . 
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The damage criteria , and to some extent the cyclic plasticity models , are novel 
and without sound experimental supporting evidence . However , it is believed that in 
combination with engineering judgment, they can be used to obtain useful qualitative 
results . 

The present in-core computer program is limited to small structural components. 
Provision for out-of-core computations would permit much broader application . 



APPENDIX-PROGRAM LISTINGS 


Following is a listing of the program CYCLIC for cyclic plasticity and fatigue 
analysis . The program includes the modifications to the NONSAP computer program 
(ref. 1) and the new overlay (number 3.8). 

Explanatory titles and descriptions of the variables used are incorporated within 
the listing. 
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3 

A 

5 

6 
7 
R 
9 

10 

11 

12 

11 

1A 

lft 

16 

17 

lft 

19 

20 
71 
2? 
23 
2A 
75 
26 
27 
7ft 

29 

30 

31 

32 
3’ 
3A 
75 
3ft 
77 

38 

39 
AO 
A1 
A 7 
A7 
AA 
Aft 
Aft 
A 7 
Ao 
A9 
c 0 
51 
* •> 

53 

S 

ftft 


* T 

♦I 

*0 

r 

r 

r 

C 

r 

C 

c 

c 

r, 

r 

c 

r 

C 

c 

c 

c 

c 

c 


OENT CYC 
NONSAP. 3 

3TAPE12#! APE13#TAPE1A» 
NUNS AP • 22 


CYCLIC 


CYCLIC PLASTICITY ANO FATIGUE ANALYSIS PROGRAM 

I. KALEV 
SEPTEMBER# 1980 

THE NONSAP PROGRAM HAS BEEN MODIFIED TO INCLUDE 

1. CYCLIC PLASTICITY MODELS# AODEO AS OVERLAY 3.8# 
(♦DECK CYCLIC)# MATERIAL MODEL 9, 

FOR 2-0 FINITE ELEMENTS (PLANE STRFSS# 

PLANE STRAIN. AXISYMMETRY) 

2. SORTED OUTPUT DATA# AS FOLLOWS# 

TAPE6-0UTPUT INCLIDES DEFLECliONS# STRESSES 
FOR MATERIAL MODELS 3 TO 8# AND THF INPUT DATA 
TAP El 2 INCLUDES PAREMETERS FOR FATIGUE ANALYSIS 
EMPLOYING MATERIAL MODELS 1#2»9 ANO 
NUMERICAL STABILITY CHECKS 
TAPE13 INCLUDES STRESSES FOR MATERIAL MODELS 1»2»9 
TAPE1A INCLUDES STRAINS AND OTHER COMPUTED RESULTS 
FOR MATERIAL MODEL 9 


r 
C 

c 
0 

r *+ ** ** ** ** ** *+ ** ** ** *♦ 

t* 

*0 N JNSAP • 67 

C . M TOT * 12000 ANO NUMEST-POOO SPECIFY THE MAX. D.O.F. AND 

r. . NONUNEaP ELEMENTS FOR THE COMPUTER CORE OF 130K 

MT3T-12000 
*0 NONSAP. 8« 

200 NUM c ST«BOOO 
r 

*T 
r 


** 


** 


56 

5051 F 

57 

1, 

5ft 

2, 

59 

5052 F 

60 

1, 

61 

2# 

62 

5053 F 

67 

1, 

6 A 

2, 

ft 5 

50 5 A F 

6ft 

1 

67 

*1 NUNS 

68 

C IN 

69 

C FXC 

70 

C I 

71 

*1 NONS 


NONSAO. A60 

OUTPUT 
WR IT E ( 1 2. 
WRITE( 12# 
WRTTF( 12. 
7023 FORMAT!// 

2021 FORMaKBX 
1# 2HFT# 6X, 
2. IX. 3HJRE 

2022 FORMAT ( 2X 

1 » 2H — # 6X. 
2# IX. 3H 

WRITE (13 
WRITE ( 1 A 
WRITE ( 1 A 
WRITE (1A 
WRITE (1A 
WR1TE(1A» 

1 AT (2 


9H 


(2 


(2 


DFPSP 


DATA AND FORMATS 
2023) KSTFP.TIME 
2021 ) 

2022 ) 

/// 10HTIME STEP .I5.20X.8HAT TIME .EIO.A/) 

. 3HNEL. 1X.3HIPT. 1 X. 2HL 0# 1 X# AHI P EL # 9X# AHDEPC# 3X# 5HSME AN. 6X 

2HSX,6X#2HSY,2X»5HALPHA»3X#3HDWE#7X#2HHP#8X#2HWP 

>6X,3HWP2.6X.3HDEE) 

. 3H . 1 X , 3H , 1X.2H— # IX. AH ,9X» AH »3X#5H >6X 

2H — .6X.2H — » 2 X , 5H .3X.3H >7X,2H — .8X.2H — 

. 6 X , 3H . 6X, 3H ) 

,2023) KSTEP.TIME 
,2023) KSTEP.TIME 
,5051) 

,5052) 

# 5053) 

505 A ) 

X.3HNEL. 1X.3HTPT ,1X» AH L0.1X.6H RATIC#1X#9H YLD»1X 

WP1,3X,9HSTRA1N( 1)#3X,9HSTR A IN ( 2 ) # 3X# 9HSTP A AN ( 3 ) , 3 X 
(A)#1X#6HAL1(1)»1X# 6HAL 1 (2) , 1 X# 6HAL 1(3)»1X#6HAL1(A)) 
X,3HNEL,1X,3HIPT,1X, AH1EEL.1X.6H C0EF.1X.9H YMAX#1X 

DWP, 3X#9H0ELEPS( 1) , 3 X, 9HD EL EP S ( 2 ) » 3X# 9HD ELEPS ( 3)#3X 
(A)#1X,6HAL2(1)#1X#6HA12(2)#1X,6HAL2(3)#1X#6HAL2(A)) 
X.3HNEL.IX.3H1PT.1X. AH M.1X.6H WE,1X,9H WE2#1X 

DF » 3 X , 9H DEPSP ( 1 ) *3X» 9H DEPS P ( 2 ) » 3X# 9H DEPSP(3)»3X 
(A)#1X,6HAL3(1). 1X#6HAL3(2)#1X#6HAL3(3)#1X,6HAL3(A) ) 

, 115H ~~ 


CASE PRINT OUT UF 0FLECT10NS TU BE OMITTED FOR ALL LOADING STEPS 
EPT NU. 1 AND NO. 20, FOR EXAMPLE, MAKE THE FOLLOWING EFFECTIVE 
F ( KSTFP.NE .1.0P.KSTEP.NE.2U) GOTO 1001 
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(13*2022) 
( 13* 2024 ) 
(13*2026) 


PRINTED ON TAPE13 


PRINTED ON TAPE13 


1001 CONTINUE 
*D NONSAP. 528 

2020 FORMAT (46H PRINT OUT FCR TINE S 
*0 TDFE.290»TDFE.294 

C . STRESSES OF MODELS 1 AND 2 ARE PRINTED ON TAPE13 

WRITE (13*2020) NG 
IF ( ITYP20.E0.0) WRITE (13*2022) 

IF ( ITYP2D.E0.1) WRITE (13*2024) 

IF ( ITYP2D.E0.2 ) WRITE (13*2026) 

WRITE (13*2030) 

*0 T0FE.319 

C . STRESSES OF NUOELS 1 ANO 2 ARE PRINTED ON TAPE13 

WRITE (13*2035) N 
*0 TDFE.360 

r . STRESSES OF MODELS 1 ANO 2 ARE PRINTED ON TAPE13 

WRITE (13*2040) I, STRESS* P1*P2* AG 
*0 TDFE.^02 

C . STRESSES OF MODELS 1 AND 2 ARE PRINTED ON TAPE13 

WRITE (13*2040) I PT* STR t S S» PI » P 2* AG 
* T TDFE.415 

IF (MODEL. EO. 9) GOTO 5040 
♦I TDFE.419 

■?040 CONTINUE 

C . HEADLINES FOR STRESSES OF MODEL 9 ARE PRINTED ON 

WRITE (13*2020) NG 
IF ( ITYP2D.E0.0) WRITE (13*2022) 

IF ( 1TYP2D.E0.1 ) WRITE (13,2024) 

IF ( ITYP20.EQ.2) WRITE (13*2026) 

*0 TDPE.491 

2020 FOR MAT ( / / //4 6H STRESS CALCULATION! 
♦0 NATRT2.74 

9 WRITE ( 6* 2501 ) ( PROP ( I ) * I- 1*NCGN ) 

WRITE (6*2061) 

RETURN 


PRINTED 


TAPE13 


RETURN 
MATRT2.137 
? 5 01 FORMAT ( 1H 


(1H »4X,*2HE 
* 1H , 4 X* 4 2H VNU 
*1H *4X* 42HYT1 MISES 
* 1H * 4X* 4 2HE T 1 SLOPE 
* 1H , 4X* 42HYT2 MISES 
»1H *4X*42HET2 SLOPE 
, 1H , 4X* 42HYT3 MISES 
*1H »4X* 42HET3 SLOPE 
,1H »4X,42HYC1 MISES 
* 1H ,4X, 42HEC1 SLOPE 
* 1 H * 4X, 42HYC2 MiSES 
* 1H *4X, 42HEC2 SLOPE 


PROP ( 
PROP ( 

1ST SURFACE,! L OADING. . PR OP ( 


SLOPE 2ND 
MISES 3RD 
SLOPE 3RD 
MISES 1ST 


1ST SURFACE, 1 
2ND SURFACE, 1 
2ND SURF»CE*1 
3RD SURFACE*! 


LOADING. ,PROP( 
LOADING. .PROP( 
LOADING. .PROP( 
LOADING. .PRUP( 


_ _ SURFACE, 1 LOADING. .PPOP( 

1ST SURFACE, RELOADING. .PRDP( 9). 
1ST SURFACE, RELOADING. . PROP ( 10) ■ 
2ND SURFACE, RELOADING. .PROP(ll)' 
2ND SURFACE, RELOADING. .PR0P(12)- 
3RD SURFACE, RELOADING. .PROP< 13)- 
3RD SURFACE*RELOADING..PROP(14). 


) * 1H 

: , 1H 

F OR MAT ( 


, 1H *4X, A2HYC3 


SLOPE 

MISES 


> 4 Y, 4 2HEC 3 SLOPE 


»4X*42HRULE 15 ) > 


,4X,45H IF RULE-0. RIGID PLASTIC /* 
*4K* 45H IF RULE-1. ISOTROPIC HARDENING /* 
, 4 X, 45H IF RULE-2.00 PRAGER KINEMATIC HARDENING /, 
* 4X, A 5H IF RULE-3.00 ZIEGLER KINEMATIC HARDENING/, 
,4X,45H IF RULE-4.00 MROZ KINEMATIC HAROENING/* 
,4X, 45HC0MBINED R UL E- . X X ( ISOTROPIC ) + ( 1-. XX ) K I NEMATI C ) 


RIGID PLASTIC /* 
ISOTROPIC HAROENING /* 
PRAGER KINEMATIC HARDENING /, 
ZIEGLER KINEMATIC HARDENING/, 


MATRT2.97 

D 29H FQ.9, CYCLIC PLASTICITY 

1NITWA.63 

11 CALL OVERLAY (4HNSAP, 3* 10, 6HRFCALL) 
INITwA . 65 

12 CALL OVERLAY (4HNSAP, 3, 11, 6HRECALL ) 
STSTN.181 

11 CALL 0VERLAY(4HNSAP,3, 10, 6HRECALL) 
STSTN. 184 

12 CALL 0VERLAY(4HNSAP*3, 11, 6H RECALL) 
0VL38.2 

OVERLAY (NSAP, 3, 10) 

0VL39. 2 
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0VERLAY(NSAP>3» 11) 


14? 

144 

145 

146 

147 

148 
140 
1 50 
151 
15? 
15? 

154 

155 
155 
157 
15 n 

150 

150 

1.51 

15? 

15? 

154 

155 
165 

157 

158 
150 

170 

171 
177 
17? 

174 

175 
175 

177 

17? 
170 
1 BO 
181 
1 8 ? 
1 8 3 
184 
1 ? 5 
186 
1*7 
18° 
1 *0 
100 
101 
17? 
17? 

104 

105 

175 

107 

108 
100 
?00 
201 

? 

20 ° 

204 

205 

206 
207 
20 8 
200 
?10 
?11 
?1 2 
21 ? 


C 

*0 6LT2D0.2, ELT2D7.7 
r 


r 


♦OPCK 

C 


CYCLIC 

PROGRAM ELT2D7 

COMMON /EL/ INO, ICOUNT, NP AR { 20 ) , NUME G, NEGL , NEGNl» IMAS S, IDAMP» ISTAT 
1 »«OOF*KLIN,IEIG, IMaSSN,IOAMPN 

COMMON /DIMEL/ N101,N102» N103»N104>N105*N106*N107»N108,N109,N110* 

1 N111,N112# Nll3>N114,N12G,N121,N122>N123fN124,N125 

COMMON /MATMCO/ aTRE Si ( 4 ) , S TR A IN ( 5 )> D ( A, 4 ) , IPT, NE l 
COMMON A ( 1 ) 

0IMENS1LN TM1I 
EQUIVALENCE ( N P A R ( 1 0 ) , NI NT ) 

EQUIVALENCE (A#IA1 

FOR ADDRFSSES N101,N102,N1Q3, . . . SEE SUBROUTINE TCDMFE 
IF (TNO.NE.O) GO TO 100 

INITIALIZE UA WORKING ARRAY 


IOW « 27 

NPT» NT NT*NJN T 

NN*N 11 0 + (NEL - 1J*NPT*IDW 
MATP*Itt(N107 + NEL - 1) 

N M ■ N 1 0 9 + t MAT 0 - 1 ) + 4 

CALL ICYCLIc (A<NN)»A(NN>,A(NM),NPT> 

RETURN 


C 


C 


r 


r 

C 

r 

r 


r 

C 


FIND STRESS-STRAIN LAW AND STRESS 

100 IPW-27 

NPT*NiNT*NINT 

NN-N110 ♦ (NEL - l)*NPT*iDW ♦ (IPT - 1)*IDW 
M ATP“I A (N107 ♦ NEL - 1) 

NM* N 109 + ( M AT P - 1>*4 

CALL CYCLIC (A(NM),A(NN)» A ( NN+4 ) , A (NN+8 ) » A (NN+12 ) . A(NN+16) 

1, A ( NN + 20 ) » A ( NN+21 ) , A(NN+22>, A ( NN+23 ) > A ( NN + 24 ) , A ( NN + 25 ) . 

2 A( NN+2fc) ) 

RETURN 

END 

SUBROUTINE ICYCLIC < WA » I WA » PP GP , NPT ) 

DIMENSION WA(27,1),1WA(27,1),PR0P(1> 

SET iNITIAL STRESSES AND STRAINS TO ZERO 
SET INITIAL YIELD POINT TO PROP(?) 

DO 10 J ■ 1 » NP T 
DO 15 I *1 » 20 
15 WA( I, J ) «0.0 

WA<21, J ) ■ PROP ( 3 ) 

WA ( ?2» J ) “PROP ( 3 ) 

IWA ( 23 » J) * 0 
I WA ( 24 » J ) ■ 0 
W A ( 2 5* J ) *0 . 

WA( 26, J ) “0. 

10 WA(27.J)«0. 

RETURN 

END 

SUBROUTINE CYCLIC ( PROP. S I G, E PS » AL i, AL2, AL 3 » YI E LD, YM AX, IPE L» LO, 
1WE2, WP2,DEPC ) 
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214 

215 

216 
217 
21 4 
21<J 
220 
’21 
222 

223 

224 
?25 
226 
227 
2?<» 
220 

230 

231 

232 

233 

234 
23 5 

23 6 
23’ 
238 
230 

240 

241 

242 

243 

244 

24 5 
244 
247 
24 q 
24 9 
’50 
251 
25’ 

253 

254 
’5 5 
254 
’57 
25 * 
’50 
263 
261 
26’ 

263 
26* 
265 

264 
267 
26* 
260 

270 

271 
”2 

273 

2’4 

275 

276 

277 
27* 
2 7 0 
280 
281 
282 
’83 
2*4 


C 

c 

c 

c 

c 

r 

C 

0 

c 

c 

s 

r 

C 

C 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

r 

c 

0 

c 

r 

C 

r 

r 

r 

r 

C 

c 

c 

c 

c 

r 

C 


1ST HUMBER OF STRESS COMPONENTS 

ISR NUMBER OF STRAIN COMPONENTS 

EPS STRAINS AT THE END OF THE PREVIOUS UPDATE 

STRAIN TOTAL CURRENT STRAIN 

OELEPS-STRAIN-EPS TOTAL STRAIN INCREMENT 

DEPS ■(1-RATI0)/M*DELE PS 

DEPSP PLASTIC STRAIN INCREMENT PER M STEP 
FOP PRINTING ONLY DEFSP-TOTAL OF ALL M STEPS 
RATIC PART OF STRAIN INCREMENT TAKEN ELASTICALLY 
RATIO IS APPLIED IN THE ELASTIC-PLASTIC TRANSITION STEP ONLY 


INCREMENT IN STRESSES, ASSUMING ELASTIC BEHAVIOR 
STRESSES AT THE END OF THE PREVIOUS UPDATED STEP 

--- -- TNTING 

THE STEP , THEN UPDATED 
STRESS-TAU 


DEL S IG 

SIG 

STRESS CURRENT STRESS FOR PR 
TAU * S IG AT THE BEG1NING OF 
AT THE END OF THE STEP 
SMEAN MEAN STRESS 
M NC. OF INCREMENT INTERVALS 
C 0R ELASTIC STATE M-l, EL aSTOPLAST IC STATE M-4 
FOR TRANSITION STEP M-3 TO 15 
PROP ( 1 ) YUUNG S MLDULUS* E 
P R 0 P ( 2 ) POISSON S RATIO 

PROP ( 3 ) INITIAL YIELD STRESS IN SIMPLE TENSION 
8ROP(3)»PROP(5)»PROP(7) YIELD STRESSES IN TENSION 
PP0P(4),PRQP(6),RR0P(8) TANC-ENT MODULE IN TENSION 
PR0P(9).PR0P(11),PPUP(13) YIELO STRESSES IN COMPRESSION 
PR0P(10),PR0P( 12>,PRuP(14) TANGENT MODULE IN COMPRESSION 
PROP ( 1 5 ) ■ PUL E ,-0 RIGID PLASTIC , -1 ISOTROPIC MODEL 
*2.00 KINEMATIC, PRAGER S RULE 
-3.00 ZIEGLER S RULE 

•4.00 MPOZ S RULE 

.XX COMBINED MQOEL- 

.XX< ISOTROPIC RULE)+(1-.XX)(K1NAEATIC RULE) 
THE COMBINED MODEL IS NOT INCORPORATED IN THIS VERSION 
NPARQ7J-15, NPAR(18)-IDW-?4 

AL1»AL2,AL3 TRANSLATIONS OF THE THREE LOADING SURFACES 
AL TRANSLATION OF THE CURRENT LOADING SURFACE 
ALB TRASLATIDN OF THE LOADING SURFACE BOUNDING THE CURRENT 
SURFACE. USED FOP HROZ S RULE ONLY 
ET,YY SlGPF, YIELD STR ESS OF THE CURRENT LOADING SURFACE 
IOEL- 0 .ELATIC LOADING OR UNLOADING 

ipel-i,2,3 plastic lqadine in surfaces 1 , 2,3 

IP ECUALS TO IPEL FROM THE PREVIOUS LOADING INCREMENTAL 
STEP OR FROM THE PREVIOUS SUBINCREHENTAL STEP 
LO NUMBER OP HALF CYCLFS (REVERSALS) 

YIELO PREVIOUS *ISEi STRESS OF THE BOUNDING SUPFACE 
YLD CURRENT UPDATED MISES STRESS ,ALSU CRITERION FOR 
PLASTIC FLOW INTIATIVE 

Y H A X HTSES STRESS OF THE ROUNDING SURFACE WHEN UNLOADING 

TT IS SAVED UNTILL THE NEXT UNLOADING 
INITIALLY YLD* YT1 , IF VP.LT.O YLD-YMAX FOR ISOTROPIC MODEL 
OR YLD*( YMAX-2+YC1) FOR KINEMATIC MODEL 
WP-(TAU-Al)*DELSIG( ASSUM. ELASTIC BEHAVIOR), FOP UNLOADING 
VP1-TAU+0EPSP , WE«TAU*DEPS-WP1 , OF- ( TAU-AL ) * < TAU-S IG ) 
WP2-YLD*DFP CUMULATIVE PLASTIC WORK 
HP-YLD-YIELD/OEP , HEP-OFP/YLP 

HE2«A8S(TAU+S1G)/2*DELEPS CUMULATIVE DURING ELASTIC STAGE 

DWP- (TAU-S IG)*OE PSP , DWE»( TAU-S IG)*OELEPS-DWP 

DER EQUIVAL. PLASTIC STRAIN ,CEE EQUIVALENT TOTAL STRAIN 

.ENT PLASTIC STRAIfi INCREMENT, DEP CLMU 

.D BE ZERO 

(STIC , -1 P. STRESS, -0 P. STRAIN AND AXI 
TRAIN, -2 P. STRESS 
(STIC RELOADING. 

ELCADING WHEN VE2 .GT . . <5*R , R«2*YC1**2/E 
START OF FULLY PLASTICITY CYCLE 
ELOADING WHEN WE2.LT .0. 2*R » 

START OF FLUCTUATING CYCLE 
ELOADING WHEN W62.LE..<5*R AND .GE..2+R, 
START OF FULLY PLASTICITY CYCLE, TOO 
I PLASTIC LOADING AFTER THE FIRST 
ELOADING STEP 


n 

• 

DEPC 

CUMULATIVE EQUIVA 

c 

• 

COMP- 

0 E PS P ( 1 +2+4 ) SHOU 

r 

• 

CCEF 

FOR PERFECTLY PL 

r 

• 

ITYP2D -0 AXIS, -1 P.S 

C 

• 

IRE - 

INDICATOR FUR PL 

C 

• 

IF 

IRF-1 

PLASTIC R 

r 

• 


INDICATES 

r 

• 

IF 

IRE-2 

PLASTIC R 

C 

• 

INDICATES 

c 

• 

IF 

IRE-3 

PLASTIC R 

r 

t 



TNOICATES 

c 

• 

IF 

IRE-0 

ELASTIC u 

c 

• 

LOAOING/R 
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385 
2 86 
2«7 

289 
?90 
29 1 
29? 
?o 3 
79 4 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 
30" 

306 

307 
303 

309 

310 

311 

312 

313 

314 

315 

316 

317 
31" 

319 

320 
3?1 

322 

323 

324 

325 

326 

327 

328 

329 

330 

331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 
34? 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 


C 

c 

r 

C 

r 

r 

C 

r 

A 

C 

c 

c 

<* 


PRINTED VALUES- WP> WP 1 , WE* 0W P, DUE, DF, OE P ,DE E ARF TOTAL OF 
M INCREMENTS PER STFP 

WE2 IS CUMULATIVE FOR ELASTIC REGION ONLY 
DEPC.WP2 ARE CUMULATIVE FUR ALL STEPS 
UUTPIT DATA ARE SORTED AS FULLORS- 
TAPE6 DEFLECTIONS ( TN ADDITION TO INPUT DATA) 

TAPE12 DATA FOR FATIGUE ANALYSIS AND NUMERICAL STABILITY 

CHECK - DEPC, SMEAN»FT»SX,SY>ALPHA#DWE»HP»WP>IRE»WP2 
* DE E 

TAP E13 STRESSES 

TAPE14 STRAINS . SURFACES TRANSLATIONS » OTHER RESULTS 


CuMMQN /EL / 1ND t ICOUNT , NP AR < 20 ) » NUMEG, NE GL »NE GNL » IMA SS » IDAM P , ISTAT 
1 ,ND0F,KLIN,IEIGtIMASSN,10AMPN 

COMMON /VAR / NG, KPR l, M CD EX, K$TEP» IT E» 1TEMAX , IR EF» I EOR EF» I NbC MD 
COMMON /MATMOD/ STR ESS ( 4 ) >STR A IN ( 4 I , C ( 4» 4 ) . IPT * NEL 
COMMON /UlSDER/ DTSD(5) 

DIMENSION PROPm,STG(l),EPS(l) 

DIMENSION TAU(4),DELS1G(4)»DELEPS(4),0EPS(4),STATE(4) 

D1MFNS1CN AL2(1)»AL3(1)»DEPSP(4)»AL1(1)»AL(4)»AL8(4) 

DIMENSION CC(4,4),CP(4,4> 

EOUIVALENCF (NPAR(3),IN0NL) »( NPAR(5),ITYP2D) 

DATA SGLAST/1000/ » STATE/ 1HE* IHPjIHP.IHP/ 
wP-HP-n. 

DEPSP(l>-DEPSP(2)-DFPSP<3)-DEPSP<4)-0. 

hPI-WE-DF-OWE-DWP-COMP-DEP-O. 

IRE-0 

M-l 

RATIC-O. 

SFAC-COEF-1. 

ALI1)-AL(2)-ALC3)-AL(4)«0. 

DO 101 1*1.4 
DO 101 J-1.4 
101 CC(1,J>»0. 

IF (IPT.NF.l) GO TO 110 
YT1 • PROP ( 3 ) 

YT2-PRGPC 5) 

YT3- FRPP( 7) 

YC1-FR0P( 9) 

YC2-PRCP( 11) 

YC3-PP0P(13) 

ET1-PRUPI 4) 

FT2- PROP ( 6) 

ET3- PRO p ( 8 ) 

EC1-PRUP( 10) 

C C2-PRCP(12) 

EC3* PROP ( 14 ) 

RULE-PROP (15) 

ET-ET1 
Y Y= YT1 

IST-4 

IF (ITYP2D.FQ.2) IST-3 
ISP-3 

IF ( IT YP2D.FQ.0) TSR-4 
YM-PROP(l) 

PV-PRQP(2) 

01-PV/(PV - 1. ) 

A2-YM/C1.+PV) 

B2-( l.-PV>/( 1.-2.0PV) 

C 2 - P V / ( 1.-2.7PV) 

C1-A2/2. 

BM-YM/(1. - 2.+PVJ/3. 

IF (ITYP20.E0.2) GO TO 105 
PLANE STRAIN / AXISYMMETRIC 
B1-A2+C2 
AI-Bl+A? 

GO TO 110 

PLANE stress 

105 A1-YM/(1.-PV*PV) 
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356 

357 

358 
35<J 

360 

361 

362 

363 
3*4 

365 

366 

367 
36« 
369 
379 

371 

372 

373 

374 

375 

376 

377 

378 

379 
300 
381 

982 

383 

384 

385 

386 

387 
389 

389 

390 

391 
•jg 9 

393 

394 

39 5 

396 

397 

398 

399 

400 

401 

402 

403 

404 

40 5 

406 

407 

408 

409 

410 

411 

412 

413 

414 

415 

416 

417 

418 

419 

420 

421 
472 
429 

424 

425 
476 


C 

c 


Bl-Al+PV 

110 YLO - YIELD 

CALCULATE INCREMENTAL STRAINS 
DO 120 I* 1* I SR 

120 DELEPSm - STRAIN(I) - EPS(I> 

IF (1TYP20.E0.2) D EL EPS ( 4 1 -Dl* ( DE L EP S ( 1 ) +DELE PS ( 2 ) ) 
f AU ( 4 ) -0 • 

DO 162 I-1,TST 

162 TAUm«5IG<I > 

IF ( ITYP2D.E0.2) STR AIN ( 4 ) - EPS ( 4 ) 


OELSIG(l) - 
DELS 1G ( 2 ) • 
0 EL S IG ( 3 ) - 
DELS IG ( 4 ) - 

IF (ITYP2D. 
DELS IG ( 4 ) . 

IF (TTYP20. 
DELS IG ( 1 ) - 
DELS IG ( 2 ) - 
OELSIG ( 4 ) - 
150 TAU ( 4 ) - 0. 
IF ( I p 6 L • GE 
IF (iPEL.GE 


4 Bl+DELEPS (2 ) 
4 A1*DELEP$ ( 2 ) 


A1*0ElEPS<1) 

B1*0ELEPS(1! 

C 1*0EL F PS ( 3 ) 

0 . 

E0.2) GO TO 150 
81 * (OFLFPS( 1J+DELFPSC2) ) 
EO.l) GO TO 150 
OELSIG ( 1 1 4 B1*0ELEPS<4) 
DELS I G< 2 ) 4 B 1* DE L E PS (4) 

OE L S IG ( 4 ) 4 Al*0ELEPS<4) 

.1) "«4 

.1) GO TO 163 


TF MATERIAL IN THE PLASTIC RANGE SK A P THE FOLLOWING 
IF MATERIAL IN THE ELASTIC RANGE CALCULATE STRESSES 

DO 160 I • 1 » I ST 

160 TAU(i) « SIG<I) 4 UELSIC-m 

CHECK WHETHER *TAU* STATE OF STRESS FALLS 
OUTSIDE THE LOADING SURFACE 
VE-OWE-O. 

DO 164 1*1,4 
WE»WE4DELFPS II) ♦TAU ( I ) 

164 0WE-DWE4CELEPS (TMOELSIGII) 

DU 203 1-1,4 

703 WE7*WF2+0ELFPSCI)*ABS(TAU(I)+SIG(I))/2. 

SM»(TAUfl)+TAU(2)+TAU(4) — Al 1< 1 ) -AL1 ( 2 ) -AL1 ( 4 ) ) /3 . 

5X-TAU(l)-ALl<l)-SM 

SY-TAU(2)-4L1(2)-SM 

S5-TAU C 3 )-*L 1 < 3 ) 

SZ»TAU(4)-4L1(4)-SM 

FT 1*1. 5*< SX*SX4SY*SY4SZ*S ?42.*SS*SS) 

FT-FT1-YLD442 

WP-SX*0ELSIG(1)4SY*DELSIG(2)4SS*DELSIG(3)*2.4SZ*0FLSIG(4) 

IF (L3.GE.1.AND.RULF.GE.2.) GOTO 167 
IF (FT) 170,170,300 
167 CONTINUE 

CHECK FOR PLASTICITY RELOADING 
AVU1NTING EARLY NUMERICALLY RELOADING 
FOR FULLY CYCLIC RELOADING IRE-1 OR IRE-3 
FOP RELOADING AT THF SAME STRESS POINT IRF-2 

IF (WP.LT.O.) GOTO 170 

TF ( (FT1-YC1**2) .LE.O. ) GOTO 170 

R-2.*(YC1**2)/YM 

IF ( ABS(WE2) .GT.0.94R) GOTO 190 

IF (ABS(WE2) ,LT.0.2*R) GOTO 191 

IPE- 3 

GOTO 300 

190 IRE-1 
GOTO 300 

191 IRE-2 
GOTO 300 

STATE OF STRESS WITHIN LOADING iURFACE - ELASTIC BEHAVIOR 
170 I PEL *0 
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STRESS(4) * o. 

DO 180 I- 1 » I ST 
180 STRESS(I) - TAU(I) 

IF ( ITYF2D.EG.2) STRAIN ( 4 ) * EPS ( 4 ) ♦ Dl* ( OFLEPS ( 1 ) ♦ DFLEPS(2)) 
DO 4 60 I - 1 » I SR 
00 4 60 J » 1 > I SR 
460 C (I. J) -0. 

C(1*1)«AI 
C (2* 1) -81 
C < 1 » 2 ) - B1 
C(2»2)"A1 
C ( 3* 3) -Cl 

IF ( 1TYP2D.E0.1 ) GOTO 400 
IF ( ITYP20.EQ.2) GOTO 470 
C(1*4)-B1 
C(2*4)*B1 
C ( 4 » 1 ) ■ B1 
C (4»2)-Pl 
C (4, 4) «A1 
GOTO 400 
470 C(4»1>*R2 
C ( 4* 2) -8? 

C(4, 3) -0. 

C ( 4 » 4 ) * A 2 
GO TO 400 

STaTF up STRFSS outside LOADING SURFACE - PLASTIC BEHAVIOR 
DETERMINE FAR I OF STRAIN TAKEN FlASTiCLY 

300 IF UPEl.FO.O.ANO.IPF.NE.2) LO-LO+1 
WE2-0. 

WE-DWF.O. 

jM-(SIG(1)+S 1G(? J*$IG(4)-aL 1( l)-ALl(2)-AtlC4») /3. 

SX-SIG ( D-SM-AL 1( U 
SY-SIG(2)-S*-AL1<2) 
jS-S IG ( 3 ) -AL 1 ( 3 ) 

S7-SIG(4)->M-AL1(4) 

DM = (DEL3IG(l)+DEL$IG<2)+DEL$IG(4))/3. 

DX - DELSIGt 1 ) - DM 
OY - 0 E LS I G ( 2 ) - DM 
OS - DR LSI G ( 3 > 

DZ - DELS1G( 4) - DM 

A » OX*nx + 0 Y*OY + 2.*DS*0S + DZ*OZ 
B = SX*DX + SY*DY + 2.*SS*DS ♦ SZ+OZ 
IF (10.LT.1) ALD-YT1 
IF (LO.GE.l) ALD-YC1 

E - SX*SX + 5 Y*S Y + 2.*SS*SS + Sl*SZ - 2.*ALD*AL0/3. 

RATIO-O. 

IF (TPcL.GT.O) GOTO 3Co 
TF ( <B*B-A*E) .LT.O) GOTO 306 
RATIOM-B + SORTC B*3-A*E ) )/A 
IF (PATIO. GT.l . ) RATIO-1. 

306 CONTINUE 

DO 350 I * 1 » I ST 

850 TAU (I) - SIG(I) + RAT 10* OF LS IG ( I ) 

IF ( ITVP2D.FQ.2 ) STP AIN ( 4 )« EPS ( 4 I + R ATI Q+Dl* ( OEL EPS ( 1 ) 

1 + DEL F P S ( 2 ) ) 

IF (RATIO. EQ. 1. > GOTO 170 

DETERMINE NUMBER OF SU BI NC P EMENTS- M. AND STRAIN INTERVAL 
IF ( FT.LT.O) GOTO 307 
M-20.*SQRT(FT)/YL0 + I 

307 CONTINUF 

IF (M.LT.3) M-3 
IF (M.GT.15) M - 1 5 
163 CONTINUE 

XM - (1. - R AT I 0 ) / M 
00 3 BO 1 - 1 r 4 

380 OEPS(T) • XM*DELEPS ( I ) 


... CALCULATIUN OF ELASTQPLASTIC STRE5SES (START) 

LOOP FuR M SUBINCREMENTS. AT THE FIRST LOOP YLO- YIELD* 
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498 

C 

• 

TaU-SIG, and the computed plastic stiffness is identical 

• 

499 

C 

• 

TO THE ONE AT THE END OF THE PREVIOUS LOADING INCREMENTAL 

• 

*!00 

C 

• 

STEP 

• 

501 

c 




502 



DO 600 IM-l.M 


503 

C 

• 

IF UNLOADING , SAVE COMPUTATIONS 

• 

504 



IF (WP.LT.O.AND.IM.GT.l) GOTH 600 


50 5 



IF <PROP(4).EQ.O.) GOTO 041 


506 



IP- IPEL 


507 



IF (IPEL.LT.l) IP-1 


50« 



IPEL-3 


500 

c 

• 

FIND THE CURRENT LOADING SURFACE WHEN LO-1 

• 

510 



IF (LG.GT.l) GOTO 993 


511 



TF (YLD.LT.YT3) IPEL-2 


512 



IF (YL0.LT.YT2J IPFL-I 


513 

c 

• 

IPEL CANNOT DECREASE DURING LOADING 

• 

514 



IF (IPEL. LT. IP) I PEL* I P 


515 

c 




516 



IF (IPEL-2) 901*902,903 


517 


001 

ET-ET1 


513 



YY* YT1 


510 



GOTO 994 


520 


002 

ET-ET? 


521 



YY-YT2 


522 



GOTO 994 


523 


003 

ET- ET3 


524 



YY-YT3 


525 



GOTO 904 


526 


003 

CONTINUE 


527 

f* 

• 

FIND THE CURREN1 LOADING SURFACE WHEN LO ,GT.l 

• 

523 



IF(TRF.E0.2.AN0.LU.FQ.l) GOTO 964 


520 



IF( IRE .60.2. AND .LO.GT • 1 ) GOTO 065 


530 



IF ( YLD.LT.ABS( YMAX-2*YC3 ) ) IPEL-2 


531 



IF ( YLO.LT .ABS( YMAX-2*YC2 ) ) iPEL-l 


532 



GOTO 966 


533 


064 

IF (YLD.LI.YT3) IPfL-2 


534 



IF (YLD.LT.YT2) IPEL-l 


535 



GOTO 966 


536 


065 

IF (YLO.LT. YC3) IPH-2 


53^ 



IF (YLO.LT. YL2) IPEL-l 


538 


066 

CONTINUE 


530 



IF { RULE . GT. 1 ) GOTO 915 


540 



IF (YL0.Lf.ABS(YMAX+2.*YC2-2.*YCl)) JPEl-1 


541 



IF (YLD.LT.ABS(YHAX+2.*YC3-2.*YC1)) IPEL-2 


542 


015 

CONTINUE 


543 

r 

• 

IPEL CANNOT DECREASE DURING LOADING 

• 

544 



IF (IPEL. LT. IP) I PEL -I P 


545 

c 




546 



TF (IPEL-2) 911,912,913 


547 


oil 

ET-FC1 


543 



YY-YC1 


540 



GOTO 994 


550 


012 

ET-EC2 


551 



YY* YC2 


557 



GOTO 994 


553 


013 

ET-EC3 


554 



YY-YC3 


555 


004 

CONTINUE 


55A 

r 

« 



5*7 



D2-YM*ET/ (YM-ET) 


553 

c 




550 



SET SURFACE TRANSLATIONS 


560 

c 




561 



IF (LO.GT. 1) GOTO 5004 


*62 



G1-YT1 


563 



G2-YT2 


564 



G3-YT3 


565 



GOTO 5007 


566 

5004 

CONTINUE 


567 



G1-YC1 


568 



G2-YC2 
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G3-YC3 

5007 CONTINUE 

IF (IPEL-2) 981*982*983 
981 DO 939 J= 1 * 4 
9P9 AL ( J ) ■ 4 L 1 ( J 1 
DU 5020 1-1*4 

5020 ALB ( I ) * AL2 ( I ) 

YRA-G2/G1 
GOTO 986 

902 00 968 J- 1* 4 
0 8 0 AL ( J )-AL2( J) 

DO 5021 1-1,4 

5021 ALB( I)«AL3(I ) 

YRA-G3/G2 

C . TC INSURE TANGENCY OF THE SURFACES 

IF (IPEL.EO.IP) GOTO 986 
DU 5001 J « 1 > 4 

'iOOl AL( J >-TAU< J)-G2/G1*(TAU(J )-ALl(J)I 
GOTO 986 

983 DO 907 J-1.4 
087 AL ( J ) * A L3 ( J ) 

YPA-1. 

DU 5003 1-1*4 
F903 ALB( I)-TAU(I) 

C . TO INSURE 1ANGENCY OF THE SURFACES 

IF (IPEL.EO.IP) GOTO 986 
DO 5002 J - 1 * 4 

50 02 AL( J >«TAU( J ) -G3 /G2* < TAU ( J )-AL2( J ) ) 

986 CONTINUE 
C 

C FOP US THE ELASTO-PLASTIC MATEP1AL HATP1X 

C 

HPRIMF-2.+D2/3. 

BETA-1. 5 /YY/YY/U.+HPRIME/AZ) 

IF ( RUL f . E Q • 1 ) BETA-BETA*YY*YY/YLD/YLD 

941 IF ( RULE • EO . 0. ) BETA-1. 5/ YT1/ YT1 
BETA1-BFTA 


IF (RULF.GE.2) GCTO 305 
DU 715 1-1*4 
715 AL ( I ) - 0 . 

305 CONTINUE 

SM-( (TAU(1)-AL(1))+(TAU(2)-AL(2))+(TAU(4)-AL(4) ))/3. 
SX»TAU(1)-AL(1)— SM 
SY-TAU(2)-AL (2)-SH 
SS-T AU( 3 )-AL (3) 

SZ«TAU(4)-AL (41-SM 

CHECK FOk UNLOADING IF WP.LT.O. 

WP-SX*OELSIG ( 1 ) +SY*DELSIG (2 ) + SS*DELS IG( 3 ) + 2 • +S Z+DELSIG( 4) 
IF (WP.LT.O.) BETA-0. 


C(l» 1) » A2 * (B2 

C ( 1* 2 ) - A2 ♦ ( C2 

C ( 2* 1)-C<1»?) 

C ( 1* 3) - A2 * ( 


C ( 3* 1 ) -C (1* 3 ) 
C (2.2) - A2 


A2 * (B2 

C ( 2 * 3 ) - A 2 * ( - BETA*SY*SS ) 

C(3*2)-C(2*3) 

C ( 3, 3 ) - A2 ♦ (.5 - BETA + SS+SS ) 

C ( 4 * 1 ) - A2 * (C2 - 8ETA*SX*SZ) 

C ( 4 * 2 ) - A2 * (C2 - BETA*SY*SZ ) 

C ( 4* 3 ) - A2 * { - BETA*SZ*SS) 

IF (ITYP2D.E0.1) GOTO 5030 
C(1*4)«C(4,1) 

C ( 2* 4) -C (4*2 ) 

C ( 3 » 4 ) «C ( 4* 3 ) 

C ( 4* 4 ) - A2 * (B2 - BETA*SZ*SZ ) 


- RETA*SX*SX) 

- BETA X* S Y ) 

- BETA*SX*SS) 

- RETA+SY+SY ) 

- BETA*SY*SS> 


IE ( ITYP2D.E0.0) GOTO 5030 


23 



640 

641 
64? 

643 

644 

645 

646 

647 
64 3 

649 

650 

651 
65? 

653 

654 

655 

656 

657 
65* 

659 

660 
661 
667 

663 

664 

665 

666 
667 
66 * 

669 

670 

671 
67? 

673 

674 

675 

676 

677 
676 
670 
680 
6*1 
68 ? 

683 

684 

685 

686 
607 
6*0 
689 
6^0 
691 
69? 
697 

696 
695 
60 6 

697 

698 

699 

700 

701 

702 

703 
706 
70 5 

706 

707 
70 8 

709 

710 


I PLANE STRESS / MODIFY OP MATRIX 

00 717 1*1,3 
A»C< I»4)/C(4,4) 

DO 717 J*I>3 

C(I»J)*C(I»J) - C(4»J)*< 

717 C(J,1> - C(I»J) 

DEPS(4)-{-C(4,l)*0EPS(l)-C(4,2)*0EPS(2)-C(4,3)*DEPS(3))/C{4,4) 
IF (WP.LT.O.) DEPS(4)-D1*(DEPS(1) ♦ DEPS(2») 
STRAIN(4)*3TRAINC4) + 0EPS(4) 

5030 CONTINUE 


C 

r 

C 


561 


* 60 
1 93 


711 

200 

123 

714 


732 

124 

773 

5012 

734 


CALCULATE ELASTIC-PLASTIC STRESSES 

IF (WP.LT.O.) GOTO 193 
DO 561 7*1,4 
CC ( I » 2 ) *0 . 0 
DO 560 1*1, 1ST 
00 560 J-l.TSR 

CC (1,2 )*CC (I»2)+C(I*J )*DEPS (J ) 

TAU(I) * TAU(I) + C(I,J) * DEPS(J) 

CONTINUE 

CALCULATE PLASTIC STBAIN INCREMENT 

IF (WP.LT.O.) 8ETA-BETA1 
C P ( 1> 1 ) «0ETA*SX* sx 
CP(1,2)»0ETA*SX*SY 
CP(1>3)-8FTA*SX*SS 
C P ( 1 » 4 ) -0ETA*SX*SZ 
CP(2»1)*CP(1»2) 

C?< 2»2 )»BFTA*SY*SY 
C D ( 2 • 3 ) *eETA*$Y*SS 
CP( 2,4 ) *8ETA*SY*SZ 
CP(3,1 )-CP(l,3) 

CP(3,2)»rP(2>3) 

CP(3»3)*8FTA*SS*SS 

CP(3,4)*PETA*S5*SZ 

C?(4»1)*CP(1»4) 

CP(4,2)«CP<2.4) 

CP(4,3)-CP(3,A) 

CP(4,4)»8ETA*SZ*SZ 
DO 711 1*1.4 
OE PS P ( I ) *0. 

DO 123 1*1,4 
DO 123 J-1,4 

DEPSP( I )*DEPSP( I ) + CP( I»J)*DEPS (J ) 

CONTINUE 

CALCULATE SURFACE TRANSLATIONS INCREMENTS 
TF ( WP .LT.O. ) GOTO 731 

A6»5X*CC.(1»2)+SY*CC(2>2)+SS*CC(3»2)*2.+SZ*CC(4,2) 

IF (RULE. IT. 2) GOTO 731 
IFUULE-3 ) 732, 733, 734 
CGNTINUP 

PRAGER HARDENING RULP 
Od 124 1*1,4 

AL ( I )■ AKI) +HPRIME*GE PSP ( I) 

AL13 >*AL(3)-HPRIME*0EFSP( 3)/2. 

GOTO 731 

ZIEGLER HAROEMNG RULE 
CONTINUE 

A12»SX*(TAU(1)-AL(1))+3Y*(TAU(2)-AL(2) ) +SS* (TAU ( 3 ) -A L ( 3 ) ) * 2 . 
1 + SZ* (T AU(4)— AL(4) ) 

A7-A6/ A12 
00 5012 1*1,4 

Aim*Al(i)+(TAU(I)-Al(I) ) * A7 

GOTO 731 

CONTINUE 

mroz hardening rule 

THE THIRD(LAST) YIELD SURFACE IS ASSUMED TO TRANSLATE 
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711 

712 

713 

714 

71 <5 

716 

717 

710 

719 

720 

721 
72? 
723 

72 4 

726 

776 

727 
770 

729 

730 

731 

73 7 
733 
736 
736 

736 

737 

738 

73 9 
76 0 
741 

74 7 
74 3 
744 
746 

746 

747 

74 0 
74*7 

750 

751 

75 2 

753 

754 

755 

756 

757 

758 
757 

760 

761 

76 7 
763 
76 4 

765 

766 

767 

76 8 

769 

770 

771 

772 

777 
774 

77 5 

776 

777 

778 

779 

780 

781 


ACCORDING TO THE ZIEGLER S Rl)LE» THUS YRA-1. AND ALR-TAU 
DO 5014 1*1*4 

5014 CC( 1,4 ) «AL3( I ) +YR A* ( TA U ( I )-AL (I) >-TAU(I) 

A9»SX*CC(1,4)+SY*CU2,4)+SS*CC(3,4)*2.+SZ*CC(4,4) 

A10-A6/A9 

AL(1)-AL<1)+A10*CC (1,4) 

AL(2)-AL(2)+A10*CC(2,4) 

AL(3)-AL(3)+A10*CC(3,4) 

AL(4)-AL(4)+A10*CC (4,4) 

731 CONTINUE 

00 712 1-1,4 

712 CC( 1,3 ) -CC < I * 3 ) +DEPSP ( I) 

CALCULATE PLASTICITY PARAMETERS 
Rl-TAU(l) 

R2- T AU ( 2 ) 

R 3-TAU ( 3 ) 

R4-TAU ( 4 ) 

PI-OERSP(l) 

P2-0EPS P ( 2 ) 

P3-DEPSP (3 ) 

P4-0EPSP14) 

COMP -COMP + OE PSP <1)4-DEPSP( 2) +DEPSP(4) 

DFP-DEP+SORT (0. 667* ( P 1+* 2 +P2** 2 ♦ P 3** 2/ 2 . +P 4+* 2 ) ) 

IF (WP.LT.O.) DEP-O. 

WT1- R1*DEPSP(1)+R2*0EPSP(2)+R3*DEPSP(3)+R4*0EPSP(4) 

IF (WP.LT.O. ) UT1-0. 

WE- WF+ R1*0FPS(1)+R2*0£PS(2)+R3*DEPS(3)+R4*DEPS(4)-WT1 
WPl-WPl+WTl 

I UPDATE SURFACE TRANSLATIONS 

IF (HP.LT.O) GOTO 904 

IF (FRuP(41.EQ.O.) GOTO 9C4 

A-YT1/YT2 

8- YT 1/ YT3 

F-YT2/YT3 

IF (LO.LF.l) GUTO 290 
A-YC1/YC2 
B-YC1/YC3 
E-YC2/YC3 
?90 CONTINUE 

IF (IPEL-2) 971,972,973 

971 00 979 J-1,4 
979 AL1 ( J)-AL( J ) 

GOTO 976 

972 00 978 J-1,4 
AL2( J)-AL(J) 

: . TC INSURE TANGFNCY OF THE SURFACFS 

AL1( J)=TAU(J )-A*(TAU( J)-AL2< J)) 

978 CONTINUE 
GOTO 976 

973 DO 977 J-1,4 
AL3< J)-AL(J) 

: . 10 INSURE TANGENtY OF THE SURFACFS 

AL2(J)-TAU(J)-E*(TAU(J)-AL3(J)) 

AL1 ( J ) = T AU( J )-8*(TAU(J)-Al3(J ) ) 

977 CONTINUE 
976 CONTINUE 
9 0 4 CONTINUE 

i \ UNLOADING 

IF (WP.GE.O. ) GOTO 920 
DM-(TAU(l)+TAU(2)+TAU(4H/3, 

DX* T AU ( 1 )— DM 
DY-TAU(2 )-0M 
OS-TAU ( 3 ) 

DZ* T AU ( 4 ) -DP 

YMaX-SORTJ l. 5*( DX*DX + DY*DY+2*DS*DS+DZ*0Z) ) 

IPEL-0 

THE ELASTIC STRAINS ARE ADJUSTED 
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782 

783 

784 
788 
78 6 

787 

788 

789 

790 

791 

792 

793 
7«4 
795 
798 
79? 

798 

799 

800 
801 
802 

803 

804 
808 
806 

807 

808 

809 

810 
811 
812 

813 

814 
818 
816 
"17 
818 

819 

820 
8? 1 
82’ 
A?3 
874 

825 

826 
82’ 
8? 8 
890 

830 

831 
032 

833 

834 
"3 8 
8’4 

837 

838 
830 

840 

841 

847 

843 

844 
"45 
846 
84’ 

848 

849 
880 
851 
88 ’ 


r 

r 


AH-IM 

W*M 

DO 921 I-1>IST 
DO 921 J-1*|SR 

9 21 TAU(I)»TAU(! ) + C(I»J)*(DEPS(J)-DEPSP( J)/AM)*W 
OWE ■ 0. 

DO 165 I«l»4 

WE«tfE + TAU(I)*(DEPS(I)-QEPSP(I)/AM)’MV-IH) 

165 DWE-OYE+(TAU(I J-SIG(I) )* ( DE PS (I)-DEPSP < I )/AM)*W 
DO 204 I« 1 > 4 

704 WE2-WE2 + ABS(TAUm )* ( DEPS (I)-DEPSP(I)/AM)*W 
DO 962 I « 1» 4 
962 DEPSP(I)-0. 

WP1-0EP-0. 

920 CONTINUE 

DM - (TAU(1)-HAU(2)+TaU(4))/3. 

OX « TAU(l) - DM 
DY ■ T A U ( 2 ) - DM 
DS • TA li ( 3 ) 

DZ * T AU ( 4 ) - DM 

IF ( PR0PC4 1 . EO.O. ) GO TO 580 
STRAIN-HARDENING MATERIAL - UPDATE YLD 
YLO- SORT (1.5 * ( DX*UX+OY*OY+2 .*US*US+DZ*0Z) ) 
IF (UO.LT.O.) YLD-ABS( YMAX-2*YC1) 

IF (WP.LT.O. AND. RULE. EO.l) YLu-YMAX 
GO TO 600 

PERFECTLY PLASTIC MATERIAL 

080 FTA-.5MDX+DX + OY*OY ♦ DZ*DZ) + CS*DS 
FTB» (YLD*YLO >/3. 

FT-FTA - FT9 
IF (PT.EQ.O) GO TO 600 
IF ( ITYP2D.EQ.2 ) GO TG 590 

C0EF--1. + 30RT(FTB/FTA) 

IF (WP.LT.O) COE^-O. 

TAU(l) - TAU(l) + COFF*DX 
TAU ( 2 ) - T AU ( 2 ) ♦ COEF + OY 
TAUC3) - TAU13) ♦ COFF*DS 
TAU(4)»TAU(4 ) + COFF*OZ 
GO TC 600 

590 COEF-SORT(FTB/FTA) 

IF (WP.LT.O) COEF-1. 

T AU ( I) -TAU(I )"COEF 
TAU(2)«TAU(2 )*COEF 
TAU(3)-TAU(3)*C0EF 

STRAIN(4)-STRAIN(4) > (CCEF - l.)*DM/8M 


600 CONTINUE 


390 


CALCULATION OF ELA5TCPLASTIC STRESSES 


STRESS ( 4 ) * 0. 

DO 390 I - 1 * I ST 
STRESS(I) - TAU ( I ) 

FINAL STIFFNFSS MATRIX 
IF (WP.LT.O.) BETA*0. 

SM-( (TAU(l)-AL (l))+( iAU(2)-AL(2))+(TAU(4)-AL(4) ) )/3. 
SX«TAU(1)-AL (I)-SM 
5Y»TAU(2)-AL (2)-SM 
S$«TAU< 3)-AL (3) 

SZ-TAU(4)-AL (4)-SM 


C(I> 1) - 42 * 

C(l»2) ■ A2 * 

C ( 2» 1 ) “C ( 1» 2 ) 

C ( 1, 3) ■ A2 * 

C (3, 1)-C(1»3) 

C ( 2 » 2 ) « A 2 * 


(B2 - BETA*SX*SX) 
( C2 - BETA*S X*S Y ) 

( - BETA*SX*SS ) 

( B2 - PETA*5Y*SY) 


( END ) 
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BETA*SY*SS) 


053 
854 
868 
85-S 
057 
88 0 
859 
869 
061 
86? 
063 
864 
866 
866 
867 
86 8 

869 

870 

871 
07 ? 

073 

074 

875 

8 7 6 
877 
070 
879 
089 
801 
08 ’ 
089 
804 
«85 
886 
807 
000 
089 

890 

891 
89? 
093 
894 
89 c 
89* 

097 

098 

099 

900 

901 
90? 
99 3 

904 

905 

906 
90? 
90” 

909 

910 

911 
91? 

913 

914 
91 5 

916 

917 
91 3 
919 

9? 0 
921 
02? 
9? 3 


C 

c 


r 

c 


C 


r 


r 

C 


79 ? 

791 

400 

716 


166 


410 

420 


700 


C(3?2)«C<2?3) 

C ( 3# 3) - A2 * (.5 - BETA*SS*SS ) 

C ( 4* 1 ) - A 2 * <C2 - BETA*SX*SZ ) 

C(4,2) - A2 * ( C2 - RETA*SY*S7) 

C ( 4 • 3 ) * A2 * ( - 6FTA*SZ*SS) 

IF ( IT YP2D.E 0. 1 ) GOTO 791 
C ( 1? 4 ) *C ( 4*1 ) 

C<2? 4 ) »C ( 4? 2 ) 

C(3?4)-C<4,3) 

C ( 4 ? 4 } ■ A2 * (B? - BETA*SZ*SZ> 

IF { ITYP20.E0.O) GOTO 791 
PLANE STRESS / MOOIrY OP MATRIX 
00 792 1-1*3 
A*C ( I » 4 ) /C (4* 4 ) 

00 792 J-T.3 

C(I*J)-C(I*J) - C (4 ? J ) *A 
C ( J * I) * C ( I , J > 

CONTINUE 

CDNTTNUF 

CACULATE parameters 
00 716 I*l?4 
OE PS P ( I)-CC(1.3) 

S1«STRFSS( 1)-STG( 1) 

52- sT»ESS(?)-SIG<2) 

53- STRFSS(3)-SIG(3) 

S 4» STR E S5 ( 4 ) — S T G ( 4 ) 

ON- ( (STRESS (1 )-AL(l))+STRFSS(2>-AL(2» *STRESS(4)-AL(4))/3. 

0X*STRFSS(1)-Al<l)-0H 

DY*STRFSS(21-AL(2)-0M 

Dl*sTRLSs(3)— AL(3) 

0Z-STRFSS(4»-AL (41-OH 
DF»DX*S1+DY*S2+0$*S3*2.+DZ*S4 

OWP-Sl + CEPoP (1) +S2*r)EPSP( 2>+S3*DEPSP(3) + S4*DEPSP(4) 

IF (1PEL.EQ.0) GOTO 166 

DWE»S1*DELEPS <l)+32*QELEPf (?)+$3*DELEPSC3) +S4*DELEPS(4)-DWP 

CONTINUE 

P1-0FLEFS(1> 

P2»0FLFPS ( 2) 

P3-OELEPS(3) 

P4-0ELFPS( 4) 

DEE*SQRT(2./3.*<Pl**2 + ?2**2+P3**2/2.+P4**2) ) 

OMMsfREssm ♦ STRESS (2 ) + STRESS (4 )) /3. 

OX-STRESS(l) - DM 
DY*3TRFSS(2) - OM 
OS* ST P E Ss ( 3 ) 

D7»STRESS(4) - OM 

FT*SCRT<1. 5*<DX*DX+DY*DY+DZ*DZ+2.*DS*DS>> 

WP2-WP?+FT*DFP 
HE P = DE P /FT 

IF (DtP.NE.O.) HP«< FT-Y1ELD) /PEP 
DEPC-DEPCfDFP 

UPOATING STRESSES? STRAINS? YIELD 
SIG(4)«EPSC4)-0. 
on 410 1*1? 1ST 
SIG( I) - STRESS (I) 

DO 420 I* 1 * 1 SR 
EPS(I) - STRAIN (I ) 

YIELD « YLO 

IF ( ITYP20.E0.2 ) E PS ( 4 ) *S TR A IN ( 4 ) 

IF (KPRI.EO.O) GO TO 700 
IF ( IC0UNT.E0.3 ) RETURN 
RETURN 
CONTINUE 

PRINTING OF STRESSES 
IF (1NDNL.NE.2) GO TO 000 



926 
97 5 
9?6 

927 
02 8 
Q?9 

030 

031 
03 2 
033 
036 

035 

036 
03 7 
038 
030 
060 
061 
06 2 
063 
066 

065 
06 6 
067 
063 
060 

050 

051 
057 
053 
056 

055 

056 

057 
050 
050 
060 
061 
062 
063 

066 

065 

066 

067 

068 
960 
07Q 
071 
07? 
073 
076 

075 

076 

077 
970 
070 
080 
081 
98"’ 

983 


C IN TOTAL LAGRANGIAN FORMULATION, 

C CAUCHY STRESSES ARE CALCULATED ANO PRINTED 

<• 

CALL CAUCHY 

r 

o 0 0 CONTINUE 

CALL MAXNIN (STRESS, SX, SY,SM» 

SMEAN«(5TRESS(lH-STRESS(2) + STRESS(6))/3. 

C 

WRITE (12, 205 2) NEL , I PT, LO, IPEL , DEPC, SM EAN , FT, SX, SY, S M, DWE# HP , WP 
1,IRE,WP2, DEE 
r 

IF (WP. LT.O. AND. RATIO. NE. 0. ) GOTO 025 
GOTO 022 
025 CONTINUE 

WRITE (12,926) NEL, I PT, K STEP 
022 CONTINUE 
r 

IF (IPT.EQ.3) GOTO 1011 
IF (IPT.E0.9) GOTO 1011 

C . THE FOLLOWING DATA IS PRINTED ON TAPE 16 FOR THE ABOVE SPECIFIED 

C . INTEGRATION PUInTS OF EACH ELEMENT 

GOTO 1001 
1011 CONTINUE 

WRITE ( 16, 5055 ) NE L , T P I, 1 0, RA TI 0, YLD, WP1, ( STR A IN < I ) , T -1 , 6 ) 

1, ( ALUI ),T«1»6) 

WRITE ( 16,5055) NEL , 1 P T, I PE L, COEF, YM AX, DWP, ( DELEPS ( I ) , I - 1, 6 ) 

2, ( AL2( I ),T*1,6) 

WRITE (16,5055) NE L , IPT , M, WE, WE2, DF, ( DE PSP ( I ) , I-l, A ) 

3, ( AL 3 ( I ), 1-1,6) 

WRITE ( 16, 5056) 

1001 CONTINUE 
r 

IF (NG.NE .NGLAST) GL TO 802 
IF (NEL.Gf.NELAST) GO TO 806 
IF ( 1° T-l ) 810,808,810 
NGLAST « NG 
WRITE (13,2003) 

NEL4sT«NFL 
WRITE (13,2006) 

CONTINUE 

WRITE (13,2007) TPT,STATE(1PEL+1),STRESS(6),(STPESS(I),I«1,3) 

1, SX, SY, SM > FT 
R ETURN 

FORMAT ( 2X, 13, IX, 13, IX, 12, IX, 16, lx, El 2. A, IX, 

16 (F 7. 2, IX) ,F6.2,lX,F5.2,lX,F8.1,iX,f9.2,lX,I3,lX,2(E9.2,lX) ) 

FORMAT (1GX, 36HRSL0ADING AT THE SAME UNLOADING STEP, IX, 6HNFL-, 

113, IX»6H1PT*»I1,1X»5HSTEP*»I6»3X»16H REDUCE STEP SIZE) 

FORMAT (2X, 13, IX, 13, IX , 1 6, IX, F 6 . 2, IX, E9. 2 , IX, E 9 .2, 6E 1 2 . 6, 6F7 . 3 ) 

FOR M AT ( 2X, 117H 


302 

300 

306 

810 


NEL 


7052 
026 

5055 

5056 

7003 FORMAT ( 102H ELEMENT STPESa ’s TP ESS- XX STPESS-YY STR 

lESi-ZZ 3TRESS-YZ MAX STRESS MIN S TR ESS , 11 X , 5HYIELD / 

2 109H NUM/IPT STATE 

3 ANGLE, IX, 8HFUNCTI0N / ) 
7006 FORMAT (16/) 

2007 FORMAT (5X»I2»2X, A1,6HLASTIC, IX , 6E16 . 6, 1 X, 2E16.6, 1X,F6.2,1X, F8.2) 

F NO 


7/8/9 


END OF RECORD 


086 


6/7/879 


ENO-OF-FILE 
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Load 


Two reversals Overload 



Figure 1. Typical idealization of loading 
spectrum. 


True 



(a) Material idealized uniaxial stress- (b) Schematic representation of yield 

strain curve at its cyclic steady state. surfaces at initial condition and after 

translation of first surface (dotted line) . 

Figure 2. Relationship between material uniaxial curve and two-dimensional stress 
field. 
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Figure 3. Incremental translations , da, representing 
three hardening rules, a ^ and a 2 represent total 

translational components of surfaces; 0 ^ and a 2 

represent stress components . 



(a) Material uniaxial relation- (h) Coffin-Manson low-cycle 

ship between stress amplitude fatigue data of material uniaxial 

and plastic strain amplitude unnotched specimen, 

at cyclic steady state. 

Figure 4. Required input data for present approach. 
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Finite element 



Figure 5. Location of discrete points in front of crack tip 
for calculation of crack growth rate. 



Figure 6. Pair of reversals count. 
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Taj Geometry of panel. 


(b) Finite element model for one- 
quarter of panel. 


Figure 7. Example of a cracked panel under uniaxial loads. 



Figure 8. Idealized material curve in uniaxial cyclic 
steady state for cracked panel. 
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Figure 10. Numerical results for cracked panel of effect of tensile 
overloads on crack growth rate . 





(b) View A. 



(c) View B . 

Figure 11. Details of an aircraft integral stiffened skin. 





Thickness Change 



(a) Finite element model. 


Nominal Stress 


21.7 kg/mm 2 -f 
(30.9ksi) 



Time 


(b) Applied compressive loading. 
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( c ) Material uniaxial curve. 

Figure 12. Idealization of stiffened skin example. 
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(a) Equal damage curves indicating number of reversals to 
crack initiation. 
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(b) Maximum von Mises stress distributions . 
Figure 13 . Results for uncracked stiffened skin. 




(a) Modified finite element models due to crack growth. 


Case 1 Case 2 



2N=i000 , 6300 , 80000 2N=3400 , 27000 

(b) Distributions of equal damage curves. 
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(c) Crack growth rate and orientation. 
Figure 14. Results for cracked stiffened skin. 
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